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Abstract 

In this paper we study the linearized inverse problem associated with imaging of reflection 
seismic data. We introduce an inverse scattering transform derived from reverse-time migration 
(RTM). In the process, the explicit evaluation of the so-called normal operator is avoided, while 
other differential and pseudodifferential operator factors are introduced. We prove that, under 
certain conditions, the transform yields a partial inverse, and support this with numerical simula- 
tions. In addition, we explain the recently discussed 'low-frequency artifacts' in RTM, which are 
naturally removed by the new method. 

1 Introduction 

In reflection seismology one places point sources and point receivers on the earth's surface. A source 
generates acoustic waves in the subsurface, which are reflected where the medium properties vary 
discontinuously. In seismic imaging, one aims to reconstruct the properties of the subsurface from the 
reflected waves that are observed at the surface [TU1 [3J IH] • There are various approaches to seismic 
imaging, each based on a different mathematical model for seismic reflection data with underlying 
assumptions. In general, seismic scattering and inverse scattering have been formulated in the form of 
a linearized inverse problem for the medium coefficient in the acoustic wave equation. The linearization 
is around a smoothly varying background, called the velocity model, which is a priori also unknown. 
However, in the inverse scattering setting considered here, we assume the background model to be 
known. The linearization defines a single scattering operator mapping the model contrast (with respect 
to the background) to the data, that consists of the restriction to the acquistion set of the scattered field. 
The adjoint of this map defines the process of imaging in general. The composition of the imaging 
operator with the single scattering operator yields the so-called normal operator, the properties of 
which play a central role in developing an inverse scattering procedure. 

There arc different types of seismic imaging methods. One can distinguish methods associated with 
the evolution of waves and data in time from those associated with the evolution in depth (or another 
principal spatial direction). The first category contains approaches known under the collective names 
of Kirchhoff migration 5 J or generalized Radon transform inversion, and reverse-time migration (RTM) 
[37l|47l[29l[Tl|4T]; the second category comprises the downward continuation approach [TT | ITO } 14} [34 l 126] 
possibly applied in curvilinear coordinates. The analysis pertaining to inverse scattering in the second 
category can be found in Stolk and De Hoop [3211101 • The subject of the present paper is an analysis 
of RTM-based inverse scattering in the first category, with a view to studying the reconstruction of 
singularities in the contrast. As was done in the analysis of Kirchhoff methods [3S] [351 |3S], we 
make use of techniques and concepts from microlocal analysis, and Fourier integral operators (FIOs); 
see, e.g., for background information on these concepts. As through an appropriate formulation 
of the wave field continuation approach, we arrive at a representation of RTM in terms of a FIO 
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associated with a canonical graph. Over the past few years, there has been a revived interest in 
reverse time migration (RTM), partly because their application has become computationally feasible. 
RTM is attractive as an imaging procedure because it avoids approximations derived from asymptotic 
expansions or one-way wave propagation. 

The study of the above mentioned normal operator takes into account the available source-receiver 
acquisition geometry. To avoid the generation of artifacts, one has to invoke the Bolker condition [21 , 
essentially ensuring that the normal operator is a pseudodifferential operator. (In reflection seismology 
this condition is sometimes referred to as traveltime injectivity condition |31j.) RTM is based on a 
common source geometry, in which case the Bolker condition requires the absence of "source caustics" , 
that is, caustics are not allowed to occur between the source and the image points under consideration 
[3"T] . We shall refer to the assumption of absence of source caustics as the source wave multipath 
exclusion (SME). Additionally, we require that there are no rays connecting the source with a receiver 
position, which we refer to as the direct source wave exclusion (DSE), and we exclude grazing rays 
that originate in the subsurface. These conditions can be satisfied by removing the corresponding part 
of the wavefield using pseudodifferential cutoffs. 

In this paper we revisit the original reverse-time imaging procedure. We do this, also, in the 
context of the integral formulation of Schneider |36j and the inverse scattering integral equation of 
Bojarski 6 . An RTM migration algorithm constists of three main parts: The modeling of source wave 
propagation in forward time, the modeling of receiver or reflected wave propagation in reverse time, 
and the applicaton of the so called imaging condition [101 13] . The imaging condition is a map that 
takes as input the source wave field and the backpropagated receiver wave field, and maps these to an 
image. The imaging condition is based on Claerbout's [5] imaging principle: Reflectors exist in those 
points in the subsurface where the source and receiver wave fields both have a large contribution at 
coincident times. 

Various imaging conditions have been developed over the past 25 years. The excitation time 
imaging condition identifies the time that the source field passes an image point, for example, using 
its maximum amplitude, and evaluates the receiver field at that time. The image can be normalized 
by dividing by the source amplitude. Alternatively, the image can be computed in the temporal 
frequency domain by dividing the receiver field by the source field and integrating over frequency, the 
ratio imaging condition. To avoid division by small values of the source field, regularization techniques 
have been applied. An alternative is the crosscorrelation imaging condition, in which the product of 
the fields is integrated over time. Later other variants have been proposed, see e.g. [27]. The 
authors of [57] use the spatial derivatives of the fields, similarly to what we find in this work. 

We introduce a parametrix for the linearized scattering problem on which RTM is based. The 
explicit evaluation of the normal operator is avoided, at the cost of introducing other pseudodiffer- 
ential operator factors in the procedure, which is, thus, different from Least-Squares migration-based 
approaches [33j . The method involves a new variant of the ratio imaging condition that involves time 
derivatives of the fields and their spatial gradients. The ratio imaging condition, albeit a new variant, 
is hence finally provided with a mathematical proof. The result is summarized in Theorem [4] As an 
intermediate result, we also obtain a new variant of the so called excitation time imaging condition in 
Theorem [5] Moreover, we also address the relation with RTM "artifacts" (Ml [HI US S3 [22 , as well 
as certain simplifications that occur when dual sensor streamer data are available. 

The seismic waves are governed by the acoustic wave equation with constant density on the spatial 
domain R" with n = 1,2,3, given by 

[ C (x)- 2 9 t 2 -A] u (x,t) = /(x,t). (1) 

Although the subsurface is represented by the half space W 1 ^ 1 x [0, oo), we carry out our analysis in 
the full space, K™. The acquisition domain is a subset of the surface K n_1 x {0}. The slowly varying 
velocity is a given smooth function c(x). The existence, uniqueness and regularity of solutions can 
be found in [55]. We use the Fourier transform: iFu(£, u>) — JJ e _l( -^' x+ "'^u(x, t) dxdt, and sometimes 
write u for Tu. 
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The outline of the paper is as follows. In section [2j solutions of the wave equation are discussed, 
starting from the WKB approximation with plane wave initial values. The (forward) scattering problem 
is analyzed in section [3j We focus on the map from the contrast (or "reflectivity") to what we refer 
to as the continued scattered field, which is the result from a perfect backpropagation of the scattered 
field from its Cauchy values at some time after the scattering has taken place. We obtain an explicit 
expression which is locally valid, and a global characterization as a Fourier integral operator. In 
section [4] we study the revert operator, which describes the backpropagation of the receiver field. 
The relation with the continued scattered field is established. The inversion, that is, parametrix 
construction, is presented in section [5] We first carry out a brief analysis of the case of a constant 
velocity. Then we introduce a novel version of the excitation time imaging condition and show that it 
yields an inversion. Following that, we present an imaging condition expressed entirely in terms of the 
source and backpropagated receiver fields, providing the RTM based linearized inversion. In section [6] 
we show some numerical tests. We end the paper with a short discussion. 

2 Asymptotic solutions of the initial value problem 

In this section, we study solutions of the wave equation with smooth coefficients. We introduce 
explicit expressions for the solution operator for wave propagation over small times. In subsection |2.1| 
we construct an approximate solution of the IVP of the homogeneous wave equation. Using the 
WKB approximation we introduce phase and amplitude functions, which are solved by the method 
of characteristics in subsections |2.2| and |2.3[ The asymptotic solution is finally written as a FIO 
in subsection |2.4| Subsection |2.5| presents the decoupling of the wave equation and general solution 
operators. Subsection |2.6| deals with the source field problem of RTM. 

2.1 WKB approximation with plane- wave initial values 

Instead of solving |TJ directly, we solve for c^ 1 u, and consider the equivalent wave equation, 

[d* - cAcKc^u) = 0. (2) 

In the later analysis it will be advantageous that cAc is a symmetric operator. We invoke the WKB 
ansatz, 

C - 1 u(x,t) = a(x,£)e iAQ(x < t} . (3) 

A straightforward calculation yields 

e -iAa^2 _ cAc ] ae iAa = X 2 a[{d t a) 2 - c 2 |Va| 2 ] 

+ i\[2(d t a)d t a + ad 2 a - 2cV(ca) • Va - c 2 aAa] (4) 
+ d 2 a - cA(co). 

An approximate solution of the form Q is obtained by requiring first that the term 0(A 2 ) vanishes, 
resulting in an eikonal equation for a, and secondly that the term O(A) also vanishes, resulting in 
a transport equation for a. We will give these equations momentarily, and comment below on the 
vanishing of terms 0(A J ) for j < 0. 

We solve ^ with plane- wave initial values: 

u(x, 0) = 0, c(x) -x 9 t u(x, 0) = e ix€ . (5) 

The role of A is here played by |£|. The WKB type solution of the initial value problem will contain 
two terms, i.e., the ansatz becomes 

c~y x, t) = a(x, t- £) c iQ ( x <*-€) + 6( x , i; £) e W x ' t; « . (6) 
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The reason is that there is a sign choice in the equation for a, leading to the eikonal equations 

d t a + c\Va\ = and d t P - c|V/3| = 0. (7) 

Here, a covers the negative frequencies and /3 the positive ones. The transport equations can be 
concisely written in terms of a 2 and b 2 . They are 

d t (a 2 d t a) - V • (o 2 c 2 Va) = and d t (b 2 d t l3) - V • (6 2 c 2 V/3) = 0. (8) 

The WKB ansatz ^ can be inserted into the initial conditions (j5j). This straightforwardly yields 
initial conditions for a, 0: 

a(x,0;£)=/3(x,0;£)=£-x. (9) 
The initial conditions for a, b can be given in the form of a matrix equation, 

/ 1 1 Wa(x,0;O\ M 

The two terms in Q are not independent. The initial value problem for a can be transformed into 
the initial value problem for /3 by replacing £ with — £ and setting /3(x, f;£) = — a(x, t; — £). Further 
analysis shows that 6(x, i; £) e 1/3 ( x '* ; £) in (K3j) is in fact the complex conjugate of a(x, t; — £) e la ( x, *'~^ . 



2.2 The phase function on characteristics 

The method of characteristics (TTJ section 3.2] will be used to solve the eikonal and transport equations, 
as usual. We first solve the initial value problem for a(y,i;£), cf. |7]) and ([9|. The same procedure 
can be applied to /3. 

The characteristic equations arc formulated in terms of (y,i), (p,w) associated with (Va,<9ta), 
and a variable g associated with a. The eikonal equation is hence given by 

F(y,t,Va,d t a,a)=0, F(y, i,p,w, ? ) = u + c(y)|p|. (10) 

The characteristic equations are then 

d M _ d M _ f-(Vc)|p|\ d9 



The only non-trivial equations are those for y and p. By (y(x, £; £), p(x, t; £)) we denote a solution 
with (y(Q),p(Q)) = (*,€). 

When a is a solution to Q, Q on some open set U C M™ +1 , and (y(-)i^(')) i s a solution to 
the first equation of (fill, where (p(-), w(-), = (V y a(y(-). *(■)). &a(y(')> *(0). «(y(0.*(0))i then 



(p(-), £*;(•), ?(•)) solve the other equations of (111, and in particular a(y(x, s; £), s; £) = a(y(x, 0; £), 0; £). 



Differentiating this identity, and using the identity (da/dy) ■ (dy/d£) = 0, which is a consequence of 



the linearization of (111, it follows that 

if y = y(x, t; £) then %a(y, t; £) = x. (12) 

To verify the local existence of solutions of (|7b , ([9|, one must derive the initial conditions for (111 
from ([T]) and ([£]) for each point y, and verify that these initial conditions are noncharacteristic, i.e. 
duiF 0. The latter is trivially the case. It follows therefore from [T7] that solutions exists up to some 
finite time locally, when 9 y x becomes singular. 

To examine the ^-dependence of the constructed solution a, we note that the initial conditions for 



(11) depend in a smooth fashion on Consequently, so does a. Furthermore, a short calculation 



shows that the function a(y,i;£) is positive homogeneous with respect to £ of degree one. 
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2.3 The amplitude function 

In this subsection, we solve for the amplitude in terms of a Jacobian of the flow of the rays. The result 
in equations ( 16 ) and ( |17[) is a manifestation of the energy conservation property. The first step is to 
carefully write equation (|8| into the form 



o t q V 



c 2 Va 
d t a 



(d t + v • V + (V ■ v))o 2 



(13) 



where we define v = — c g ^ a . We used that (dt + v • V)d t a = 0, i.e. the frequency is constant on a ray. 
The field v is associated with the rays, which satisfy 



dt 



We have = The derivative | 



(i;x)=v(y(t;x),i). 



(14) 



— det 



<9y 
<9x 



det 



tr 



is hence related to V • v as 
<9y\ -1 d dy 



dt dy 



(V ■ v) det 



dy 



(15) 



This implies that 



det(<9 x y) a 2 is constant along the ray. 



(16) 

Indeed, ( 16 ) is easily established by computing the derivative % [det(9 x y(t; x))a(y(t; x), t) 2 ] and using 
(13). From ( 16 1 it follows that a(y(t;x), t; ■) — ^ det(d x y(t; x) _1 ) a(x, 0; •). Inserting the ^-dependence 
back into the notation, and using that the map x t— > y(x, £;£) is invertible results in 

a(y,t;£) 



2c(x(y,t; 0)|^| 



det(5 y x(y,t;0)- 



(17) 



2.4 Solution operator as a FIO 

In this subsection we consider more general initial values than ([5| by considering linear combinations 
of the terms in (|6| . This results in an approximate solution operator in the form of a Fourier integral 
operator (FIO) [16, 4511351 [2D] and we will review some of its properties. Our solutions so far involve 
only the highest order WKB terms and are limited to some small but finite time. 
We consider the original wave equation with / = and the initial conditions 



*(x,0)=0, 



<9 t u(x,0) = /i 2 (x). 



(18) 



Following (JsJ) , its WKB solution for time t € I, which we will denote for the moment by SiaW^y) 
is given by a sum of two terms Si 2 {t)h 2 (y) = c(y)(S , a2 (i)/i2(y) + S\,2(t)h% (y)), with 



S a2 (*)/i 2 (y) 



(2 



1/7 e i«(y,*;0-^x^i) ft2(x)dxde 

Tr)™ J J C(X) 



(19) 



Here the subscript "a" refers to the negative frequencies, i.e. phase and amplitude functions a and a. 
Then Sb2 is defined similarly, us ing and b, and refers to positive frequencies. We recall that the 
symmetry relations of subsection 2.1 imply that S^W^ = ^2(^)^2- The construction is such that t 
can be negative. 

To argue that 5 a 2 is a FIO, we will take a closer look at its phase function, i.e., 



¥>(y, t, x, |) = a(y, *;£)-£■ x, 



(20) 



5 



and observe that it is positive homogeneous with respect to £ of degree one, as it should. The stationary 
point set is given by 

r t = {(y,x,£) g yxlxl"\{0} | x = %*(y, «;€)}• (21) 
For r f to be a closed smooth submanifold of i r xXxK™\{0}, the matrix, 

d x d£tp 

needs to have maximal rank on T t , which is obviously the case j3BJ chapter VI, (4.22)]. The stationary 
point set T t is hence a 2n-dimensional manifold with coordinates (y, £). 




The stationary point set can be understood in terms of the bicharacterstics. Definition (21 1 allows 
us to express x on T t as a function xr(y, t, £) = <9|a(y,i;£). Equation (12) implies that (y, x, £) G T t 



if and only if a bicharacteristic initiates at (x, £) and passes through (y, rf) at time t where r\ must be 
given by rj = d y a(y, t; £). If (y, x, £) G FxlxR n \{0} and t G K are such that (y, x, £) G T t then one 
has dta(y,t;£) = — c(x)|£|, since the frequency dta is constant on a ray. 

The propagation of singularities of 5* a 2 is described by its canonical relation, 

n t = {((y,r,),(x,O)6rY\0xT*X\0|x = x r (y,t,O,T7 = 9Xy.*;O}- (22) 

Clearly, II t is the image of T t under the map (y, x, £) h-> ((y, 77), (x, £)). It follows from the character- 
istic ODE that the map from (x, £) to (y, tj) is a bijection, <!> t : T*X\0 — > T*Y\0 say. The canonical 
relation is hence the graph of an invertible function. Therefore, each pair (y,£), (x, £) and (y, 77) can 
act as coordinates on T t , and on n t . We observe that <I>t depends smoothly on t. 

The effect of the FIO S&2 working on a distribution v can be explained in terms of the wave front 
set. If v G £'(X), then the wave front set WF(«) of v is a closed conic subset that describes the locations 
and directions of the singularities of v. Operator S a 2 affects a distribution v by propagating its wave 
front set by composition with the canonical relation[T6 l 124 ] l45 l [46] . From the above description of II t 
it follows that 

WF(5 a2 (t)u) C $t(WP(«)). (23) 

The pair (x, —d x ip) are referred to as the ingoing variable and covariable, and (y, d y (p) as the 
outgoing variable and covariable. The idea behind the names is that S a 2, by <i> t , carries over (x, £) of 
the ingoing wave front set into (y, rf) of the outgoing wave front set [46,, p. 334]. 

So far the highest order WKB approximation was used. The notion of symbol classes for a, b 
is needed to properly include lower order terms. By replacing a by an asymptotic sum a(x, t; £) = 
S^Lo a m-j(x, t; £), with a& homogeneous of order k in £ for |£| > 1, the error in (|gJ) can be made 
to decay as |CI _Ar f° r an y N. In other words, it becomes C°° and the approximate solution operator 
becomes a parametrix. Moreover, the exact solution operator can be written in the form of c(S a 2 + Sb2) 
by the addition to a and b of certain symbols in 5^°°, which in particular decay faster than any power 
\£\~ N (unsurprisingly, the latter additions cannot be computed with ray theory). 

Solution operators for longer times have been constructed using more general phase functions. For 
us those explicit expressions are of no interest, but we note that the FIO property, with canonical 
relation characterized by $t, remains valid, as can be seen by applying the calculus of FIO's [El 
theorem 2.4.1] to the product of several short time solution operators. 

2.5 Solution operators and decoupling 

we assumed that the functions ae 1Q and be 1 ^ propagate independently as solutions 



In subsection 



2.1 



of the wave equation. In fact, this is the result of a rather general procedure to decouple the wave 
equation [15] ■ Because the results of the decoupling will be used explicitly in section [4] we give a short 
review of it here; we will examine its relation to the solution operator S"^. 



G 



We write the wave field as the vector (ui(x, t), u 2 (x 7 t)) T = (u,d t u) T . The homogeneous wave 
equation can now be written as the following system, l st -order with respect to time. 



uA _ ( A (m 

u 2 lc(x) 2 A Oj U 



(24) 



The solution can be given as a matrix operator that maps the Cauchy data at t — 0, say (mo,i(x), mo,2(x)) t , 
to the field vector at t. 

£)- a <«te) - h ™-{™ (25 » 

Naturally it satisfies the group property S(t)S(s) = S(t + s). It is invertible by time reversal. 

To decouple the system, we define several pseudodifferential operators. Let operator B be a sym- 
metric approximation of y— c(x)Ac(x) with its approximate inverse B~ x such that B 2 +cAc, B~ 1 B—I 1 
and BB^ 1 — I are regularizing operators, i.e. pseudodifferential operators of order — oo. Although the 
square root does not necessarily have to be symmetric, being symmetric has the advantage that it 
yields a unitary solution operator, as we will see. Neglecting regularity conditions, we use symmetry 
and self-adjointness interchangeably. The principal symbols of B and B^ 1 are c(x)|£| and ^jTgi re- 
spectively. The existence of such operators is a well known result in pseudodifferential operator theory, 
see e.g. [32] . We now have the ingredients to define two matrix pseudodifferential operators A and V 

by 

1 1 \ , . fl iB' 1 



V = C ^{-iB iB) and A =^l -ifl-^sfo' ^) 

which are each others inverses modulo regularizing operators. We finally define the following two fields 
(u a (x, t), tib(x, t)) T = A(u\,U2) T ■ Note that the Cauchy data can be represented by a time evaluation 
of (w a ,Mb) T . We will use the phrase 'Cauchy data' in this way also. Omitting the regularizing error 



operators, the system (24) transforms into a decoupled system for (u a , Ub) of which the first equation, 



together with its initial value, is 

d t u a = -iBu a and w a (x, 0) = u , a (x). (27) 

By removing the minus sign it becomes the equation for u^- Let and Sb be solution operators of 
the IVPs, i.e. w a (x, t) = S a (t)uo ta {x) and similar for Sb- Therefore, modulo regularizing operators 

S(t) = v( S f ^ ) A , (28) 



which means that the original IVP (24) and the decoupled system (27) have identical solutions disre- 
garding a smooth error. Because B is self-adjoint operators 5 a and Sb are unitary, which follows from 
Stone's Theorem [13]. It can be shown that S a (t) and S h {t) with tgR are FIOs [33] . 

We turn to the relation of this matrix formalism and S\ 2 — c{S a 2 + Sb 2 ), from which we derive 
a local expression of p.s.(5 a ), the principal symbol of S a . The amplitude of S a 2 is a homogeneous 



symbol, which implies that it coincides with its principal symbol, and from its definition ( 19 ) can 
thereafter be concluded that S a 2 = P-S-(S , a A 12 ). The principal symbol of a composition is the product 
of the principal symbols of its factors [16j|46], and hence — p.s.(5 a ) 2 C (x) 2 |$| • Using the solution 



of the transport equation (17), one concludes that 



p.s.(S a )(y, x, £) = x /det(9 y x(y, t; g)). (29) 
The principal symbol of Sb follows from Sb = S a 
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2.6 The source field 

In this subsection we discuss the source problem. The unperturbed velocity is a smooth function c(x). 
The source wave is the fundamental solution of a delta function located at (x s , 0) in space-time. 

[c(x)- 2 9 t 2 - A]g(x, t) = S(x - x s )S(t) m 
<?(x,0) = 0, d t g(x,0) = Q. 

An important assumption is that 

the source wave does not exhibit multipathing (SME). (31) 

The fundamental solution can therefore be approximated by an asymptotic expansion with a single 
phase function. This can in principle be found by an application of section [2. 4| and using a change of 
phase function [TBI section 2.3]. One can show that, if |x — x s | > e for an e > and t bounded, the 
fundamental solution can be written as the Fourier integral [2] 

g(x, x S) t) = ^ J A(x, x s , u) e-^-Tfx^)) duJ ^ (32) 

with A(x,x s ,U)) G S 2 and A(x,x s ,uj) = ^]^ iji(x,x s ,w). Each term is homogeneous, i.e. one has 
Af-(x,x B , Aw) = \~^~ k Ak{x, x s ,L}) for A > 1 and |cj| > 1. This holds for n — 1,2,3. The sum means 
that for each N £ N there exists a Cn > such that 

|4(x,x s ,o;)-^J o 1 ^(x,x B , W )| < CW (1 + M)*^-". (33) 



The source is real, implying that Ak(x, x s ,(jj) = Ak(x,x s ,—uj) for all fc. In (32 1 one can also view the 
separate contributions of positive and negative frequencies. 

In part of the further analysis we will use the highest order term of the source field. There exist an 
amplitude A s (x) and a cutoff cr(u), both real and such that A (x,x s ,uj) — A B (x)a(u>) (iu;)~~9~ on the 
support of a. Function a is smooth and has value 1 except for a neighborhood of the origin where it 
is 0. We also abbreviate T s (x) = T(x, x s ). The principal term of the expansion can now be written as 

g(x,t)=A s (x)d^6(t-T s (x)). (34) 
Functions A s (x) and T s (x) will be referred to as the source wave amplitude and traveltime respec- 

"-3 ... . 11-3 

tively. Operator d t 2 denotes the pseudodiffcrential operator with symbol uj t— > (j{uj){\uj)^r . The 
approximation g(x,t) matches the exact solution in case Vc = in the limit of uj — > oo. In that case 
one would have T s (x) = S^sl and A s (x) = § , A / 87r | x c l x i for respectively n = 1,2, 3 [2|. We 



r|x— x s | ' 47r|x— x s | 

define the source wave direction vector 

n s (x) =c(x)9 x T s (x). (35) 

This vector will, for example, be used to provide insight in the microlocal interpretation of the scattering 
event. 

Source waves that arrive at the acquisition set are in the context of the inversion called direct 
waves. The negative frequency part of the wave front set of the source field is given by 

S s = {(x,t,£w) G r(lxt)\0 | (x,0 - $ t (x s ,0, £ s G K"\{0}, w = -c(x)|£|}. (36) 

It contains all bicharacteristics that go through (x s ,0) in spacetime. In the region where the Fourier 
integral (32 1 is valid, direct rays are also described by the equations t — T s (x) and ^ = |£| n s (x). The 
restriction to time t c is denoted by 

S s , tc = {(x,€) G T*X\0 | (x.tc^w) G S s }. (37) 

This will be used to describe the direct waves in the Cauchy data of the continued scattered field. 
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3 Forward scattering problem 



We consider the scattering problem and formulate the continued scattered wave field as the result 
of the scattering operator acting on the reflectivity, i.e. the medium perturbation. We start with a 
description of the scattering model, essentially a linearization of the source problem. In section 3.1 we 



derive an explicit expression for the mentioned operator. It will be used in section 3.2 to define the 
global scattering operator, of which we show in theorem [2] that it is a FIO under the conditions of the 
DSE and the SME. 

3.1 Continued scattered wave field 

Here, we introduce the scattered wave field and the continued scattered wave field. Loosely stated, 
the latter is the reverse time continuation of the former. We introduce the scattering operator that 
maps the medium perturbation to the continued scattered wave field. Theorem [T] shows that a local 
representation of the operator can be written as an oscillatory integral. 

The medium perturbation is modeled by the reflectivity function r(x). The non-smooth character 
of the perturbation gives rise to a scattered or reflected wave. We assume that 

supp(r) C D for a compact D c M™ -1 x [e, oo) and some e > 0. (38) 

The last component of x describes the depth. Because the source is at the surface, i.e. x s n = 0, the 
reflectivity is zero in a neighborhood of the source. Following the Born approximation, the scattering 



problem is obtained by linearization of the source problem (30) with (1 + r(x))c(x) as the velocity. To 



find the linearization it is advantageous to first multiply (30) with c(x) 2 . The result is 



[c\ 2 - C (x) 2 AMx,£) -r(x)2A s (x)V 1 ^-T s (x)) (3Q) 
u(x,0)=0, d t u(x,Q)=0. 



The scattered wave field it(x, t) is defined as the solution of the scattering problem ( 39 ) . We have used 
that the source wave field does not exhibit multipathing (SME) and can therefore be formulated as 
the asymptotic expansion (132]) . In the forward modeling we will use the principal term to approximate 



the source, i.e. (34). The subprincipal source terms do not contribute to the principal symbol of the 
scattering operator [33] , 

The continued scattered wave field Uh is defined as the solution of a final value problem of the 
homogeneous wave equation such that the Cauchy data at t = T\ arc identical with the Cauchy data 
of the scattered field u: 

[c\ 2 -c(x) 2 AK(x,£)=0, 

Uh(x,Ti) = u(x,Ti), cW(x,Ti) = d t u(x, Ti). 

The contributions to the scattered field entirely come to pass within the interval [To, Ti], i.e. To and Ti 
are chosen such that T s (supp(r)) c [T ,Ti]. For t > T\ one has Uh(x,t) = u(x,t) but as does and 
u does not solve the homogeneous wave equation, they differ for t < T\. We also use the decoupled 



wave fields (uh, a ,Uh,b) T = A(uh,9(Uh) T , with A defined in (26). 

The continued scattered wave field models the receiver wave field in an idealized experiment. 
Idealized here means that all scattered rays are present, even rays that do not intersect the acquisition 
set. It hence represents the scattered field by being its continuation in reverse time. The reverse time 
continued wave Geld, to be defined in section |3J models the receiver wave field. 

The scattering operator F by definition maps r to {u\ 1 ,dtU\ i ) T , and we let F a and F\> map the 
reflectivity r to the decoupled components of the continued scattered wave field «h,a an d Uh,b- To 
show that F a is a FIO we derive an explicit formulation valid for a small time interval around a 
localized scattering event. Let be a finite smooth partition on D such that X^ex/ * — 1 on D. 

Using pi as multiplication operator then 

F^) = y\ S^t-hjF^pi, (41) 
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and i*b likewise. S a is the solution operator (28). The i th local scattering event is delimited by [toi,tu], 
so T s (supp(pi)) C [toij tu]- The partition is chosen fine enough such that [toi, tu] falls within an interval 
of definition of (19), i.e. the local expression of solution operator 5 a 2. 

We write p for an arbitrary member of {pi}i<=z and [to,fi] for its delimiting interval, and derive a 
local expression of the scattering operator evaluated at t\. We will prove the following 

Theorem 1. The local scattering operator F a (ti)p can be written as an oscillatory integral. It maps 
the reflectivity r to the continued scattered wave field, that is, Uh a (y, ti) — F a (tx)pr(y) and 



Wh, a (y,*i) 



1 



ei*r(y.*i.*.e>A F (y, t!,x, i)di pr(x)dx, 



(42) 



in which the phase and amplitude function are respectively defined as 
<P T (y,*i,x,£) = a(y,h - T s (x);£) - £ ■ x, 



A F (y,i!,x,0 = (ifta(y,ti -T„(x);€))" 



, a(y,ti-r g (x);£) 
c(x) 



2A s (x). 



(43) 



Here (42) is only the contribution of pr. There is a similar statement for ith,b, which satisfies 

«h,b(y,*) = Wh, a (y,*)- 



Proof. To solve the scattering problem ( 39 1 it will be transformed into a r-parameterized family of 
IVP's. DuhamePs principle states that the solution, i.e. the scattered wave field, is given by 



u(x, t) 



«(x, i; r) dr, 



(44) 



in which for each r function u(x, i; r) is the solution the homogeneous wave equation with prescribed 
Cauchy data on t = r [171 §2.4.2]: 



[<9 f 2 - c(x) 2 A]m(x, t; t) = with t 6 R 
m(x, r; t) = 

fltfi(x, r; r) = r^A^d^ 6(t - T s (x)). 



(45) 



The continued scattered wave field is the solution of the final value problem (40). Using the 

n+l 

observation that r(x)2A s (x)9 t 2 8(t — T s (x)) = if r ^ [Tq,? 1 !], it can be found by 



Uh(x, i) := / u(x,t;T) dr with tG 



(46) 



Time integration is now over the fixed interval [To,Ti], by which Uh solves the homogeneous wave 
equation. For t>T\ the wave fields u and ith coincide. Therefore, this solves (40). 

To derive the local expression we solve the r-parameterized homogeneous IVP ( 45 1 with r replaced 
by pr and evaluate the solution at t\. Let (u a ,Ub) T = A(£t,9tu) T , then u — c(u a + u^). We apply 
solution operator S 1 ^ with initial state at time r. This gives 



l (y,t 1 ;r) = 5 a2 [pr(x)2A 8 (x)9 t 2 <5(r - T 8 (x))](y,ti - r). 



(47) 



Note that ^2 involves a relative time, i.e. the difference £i — r, which is allowed because the medium 
velocity does not change in time. Then, time is as much as absolute when it agrees with the source 
time reference. 



10 



Consider Uh, a (y,^i)j i.e. integral (46) with u replaced by u a (y,ti;r) in (47). We will eliminate r 
by integration and write the field as an oscillatory integral. With the expression ( 19 ) of S a 2 and the 
application of T s (supp(pr)) C [to>£i] one derives the following integral 



Uh,a(y,ti) 



(27t)™ 



e i Q (y,ti-r;€)-i€-x o,(y,h r; g) pr ^ 2A j^ dT^Sir - T 8 (x)) dxd|dr. 
c(x) 



"+1 "+1 

We recognize two convolutions, the integral over r and operator d t 2 , the operator d t 2 can be 

n+l 

commuted to act on e 1 "- 1 ?*^. Restricting to the highest order term, one writes d t 2 [e 1Q_1 ^' x ^] = 

(i<9 t a)T^e la_1 £' x ° , which is an application of a general result of FIO theory [151 US]- Cutoff a is 
omitted to shorten the expression. This yields 

1 



Wh,a(y,il) 



(27!-)" 



e ia -*' x (id t a)- 



c(x) 



pr{x)2A s (x) S(t- T 8 (x)) dxd£dr. 



(y,*i-T;£) 



Notation [. . . ] arg means that a, dta and a within the square brackets are evaluated in given argument. 
Explicit integration finally gives the oscillatory integral in (42), (43). □ 



3.2 Scattering operator as FIO 

Here we establish that F a (ti)p is a FIO if the direct waves are excluded (DSE). We define the global 
scattering operator irF and show that it is a FIO with an injective canonical relation, i.e. theorem [2j 
Before we proceed with the theoretical aspects of the operator, we will explain what it does. The 
stationary points of F a (ti)p are given by d^ip T — 0, i.e. d^a(y, ii — T s (x); £) — x = 0. A stationary point 
(y, x, £) has the following interpretation. The source wave front excites the reflectivity at (x, T s (x)) 
in space-time, causing a scattering event. The event emits a scattered ray from (x, T s (x)) with initial 
covariable £, which arrives at (y, ti) with covariable t] = d y <p T (y, ti,x, £). Operator F a {t\)p so describes 
the scattering event and the propagation of the scattered wave over a small distance. The distance 
will be extended by application of the solution operator, see (|41[). Using the terminology introduced 



at the end of subsection 2.4 the ingoing variable and covariable are (x, £) with £ = — d x ip T (y, ti, x, £). 
The outgoing variable and covariable are (y,T]). This means that F a (ti)p carries over (x, £) £ WF(r) 
into (y,fj) € WF(u h , a (.,ti)). 
We have 

C = -d x ^ T = d t a(y, h - T s (x); d x T s (x) + £. 

Using the source wave direction vector n s (x) = c(x)9 x T s (x) and the identity 9 t a = — c(x)|£| for the 
frequency, this yields the relation between £ and £, 

C = £ - |C| n s (x), (48) 

reflecting Snell's law. Figure [l] shows the microlocal picture of the scattering event and the scattered 



ray. Equation (48) also implies that £ • n s (x) < everywhere. This is a result of the geometry of 



the reflection event with one source. Note that (48) only holds for negative frequencies. For positive 



frequencies, i.e. considering F h , one gets = n s (x) instead. In that case £'-n s (x) > everywhere 



If (x, £) is associated with a source ray, i.e. £ = |£| n s (x), then £ = by (48). In that case there 
is no reflection. We show that away from the source rays the scattering operator F is a FIO with an 
injective canonical relation, which will be made more precise. The practical implication is that source 
wave arrivals are excluded from the data before the receiver wave field is calculated. 

The direct source wave exclusion (DSE) is the removal of the source singularities contained in S s 
from the wave front set of the continued scattered wave field. Mathematically it will be applied by 
i-families of pseudodifferential operators 7T a (t) and 7Tb(t) that act on the Cauchy data (ith,a(', t), U} 1 ^(-,t)) T . 
The symbol of 7r a (i c ) is, for some fixed t c , a smooth cutoff function on T*Y\0, being on a narrow 
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•(y.*i) 



Figure 1: Propagation of singularities at (x,T s (x)) in space-time. See equation (48). The dotted line 
represents the ray. The endpoint of the ray at (y,ii) contributes to the scattered field. Here, |£|n s 
and £ can respectively be interpreted as the wave numbers of the initial and reflected waves, and C a 
normal vector that can be associated with a reflector at x. 



conic neighborhood of 3 S) t c (cf. (37)) and 1 outside a slightly larger conic neighborhood. Furthermore, 
we assume that ir a satisfies 

7T a (f) = S a (t - t c )n a (t c )S a (t c - t), (49) 

which implies that the field 7r a ith ja still satisfies a homogeneous wave equation. The symbol 7Tb satisfies 
7Tb(*;x,£) = 7r a (i;x,-£). 

Since, in the absence of multipathing, rays define paths of shortest traveltime between two points, 
we have the following property. Let x, x g D be not identical, then 

if M"\{0} and ti > such that (x,£) = $ t .(x,£) . . 

(50) 

then |T s (x) - T„(x)| < U or (x,£) g H s , Ts (x)- 

If x and x lay on the same source ray then |T s (x) — T s (x)| = f;. 

The central result is the theorem that the composition 7r a F a is a FIO of which the canonical relation 
is the graph of an injective function. Let V s t C T*Y\0 be the zero set of 7r a (i), a conic neighborhood 
of Sa,t- With 7Tb Ft, — n a F a we present the following 

Theorem 2. Operator ir a F a defined above, is a FIO. Its canonical relation is 
A={((y,t,»/,w),(x,C)) | (y,»j)G(2T\0)\V,, il teR, w = -c(y)H, 

(x,€) = * T .(x)-t(y,»i), C = C - |€| n 8 (x) , xe^}. (51) 

The projection of A to its outgoing variables, i.e. (y,t,rj,Lj), is injective. 

We will first show that the composition Tr a (t±)F & (ti)p is a FIO. Composition 7r a F a is subsequently 



defined as the sum of local contributions, like in (41 ), and will also be called the 'scattering operator' 



The canonical relation becomes the union of the local relations. A part of the proof is put in lemma[TJ 
The operator can alternatively be defined by means of the bicharacteristics of the wave equation. The 
papers [351 131] show how this can be done, although their scattering operator does not fully coincide 
with ours. 

Proof. Because TT^tjS^t — tu) = S a (t — tii)7r a (£ij) the scattering operator can be written as 

K a (t)F a (t) = V S a {t - tiiKfaO^faOft- (52) 

Again omitting subscript i to denote an arbitrary member of I we will argue that the local scattering 
operator TT a (ti)F a (ti)p is a FIO. Then 7t a F a becomes a sum of compositions of FIOs. 



The local scattering operator is the oscillatory integral (42) in which the amplitude Ap (W3| is 



replaced by n a (ti; y, d y a) Ap. This follows from the application of pseudodifferential operator 7r a (ii), 
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its symbol denoted by 7r a (ii; •, •), on the integral [THl HB]. To be able to omit the zero set of 7T a (ti) 
from the analysis of the phase <p T we define the conic set 

W Ml = {(y,x,£) e yxJfxir\{0} I (y,r,)ey Ml , (x,€) = *T.(x)-t 1 (y,»l), xec}. 

The stationary point set of the phase function, by definition d^ip T = 0, is given by 

E tl - {(y,x,€) G (yx!xR"\{0}) \ W Ml | x = %*(y,ti -r,(x);{),i6 supp(p)}. (53) 



We observe that |9 X 9^(/3 T | = |9^3 x y T | = |c^CI by definition of £ (48). Moreover 

I^CI = |I„ - li ® n s (x)| = 1 - || • n s (x). (54) 



By the DSE, applied as the omission of W e ti in (53), the condition £ || n s (x) is never met, from 
which follows that the Jacobian \d^\ is nonsingular. This implies that the derivative d^ ytX ^d^ip T 
has maximal rank, making E tl a closed smooth 2n-dimensional submanifold. The canonical relation 
relates ingoing (co)variables (x, £) with outgoing (co)variables (y, rj) and is given by 

{((y> V), (x, C)) | (y, x, G E 4l , t] = d y ip T , c = -d x ip T } . (55) 

The relation is the graph of a diffeomorphism. We postpone the proof until after the construction of 
the global scattering operator Tv a (t)F a (t) as the local and the global arguments are basically the same. 
Therefore the local scattering operator is a FIO with a bijective canonical relation. 

The local operator will be composed with the solution operator. This gives a seamless extension 
because both operators are build on the same flow. It becomes S a (t — ii)7r a (fi)F a (ti)p, which is a FIO. 
The canonical relation is determined by the composition of relations [THl HBJ • The global scattering 
operator ir a (t)F a (t) is subsequently defined as the sum (52) of the extended local operators, of which 
the canonical relation A t is the union of the local relations ( [55] ) . 

We will argue that A t is the graph of an injection Q t : (T*D\0) \ U s -> (T*Y\0) \ V s . t that is a 
diffeomorphism onto its image. We used the zero set of n a (t) expressed in the domain of &t- 

[/ s = {(x,()efl\0 | (x,€)G7.,r.(x), C = £ - l£l n s (x) , xefl}. (56) 
The injection implies that At, for a fixed t, can be parameterized by y and r], so 

A t = {((y,7,),(x,C)) | (y,v)^(T Jf Y\0)\V s , u 

(x,€) = $r.(x)-t(y,»l), C = £ - |€| n 8 (x) , xeD}. (57) 

We now prove the existence and injectivity of 8(. Without loss of generality we assume that t denotes 
a moment after the scattering event. 

Let (x, C) G (T*D\0) \ U s be given. It can be shown that the transformation £ i— ^ £ given in (48) is 
injective on the complement of U s and thus determines a unique (x, £). By ray tracing over t — T s (x), 
i.e. mapping by * t _T s (x), one finds (y,»j). 

Let (y,rj) G (T*F\0) \ V^.t be given. This uniquely determines a bicharacteristic. By ray tracing 
backwards, i.e. by $T 3 (x)-t w hh T s (x) — t < 0, the ray goes through (x,T s (x)) in space-time. If a 
second point (x, T s (x)) is met, property (50) (SME) implies that the bicharacteristic coincides with one 
from the source. The condition (y, rj) £ V s ,t (DSE) rules out this possibility, leading to the conclusion 
that x is unique. The covariable £ uniquely follows from the ray tracing, and is mapped to £ by ( 48 ) . 
The transformation O t is therefore one-to-one. 

To prove the smoothness we analyse the scattering event around a fixed point (xo,£ ), of which 
xo G supp(p), and define To = T s (xo). Now Qt can be factorized as follows 

(x,C)^(x,0^^(x,£)^My,r,). 
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The Jacobian of t becomes the product of three Jacobians, namely 



diyvn) 

d(x,C) 



diyvn) 



d(x,t) 0(x,g) 
9(x,0 5(x,C) 



(58) 



The leftmost factor in the right hand side is nonsingular because <&i_ To is a diffeomorphism. The 
rightmost factor in the right hand side is nonsingular because the map £ i— > £ has a positive Jacobian 
(54). The transformation $ ro _T s (x) is the least obvious one. We will show in lemma [I] that it is a 
smooth bijection. Therefore <d t is a diffeomorphism onto its image. 

So far t was held fixed to simplify the presentation. Time dependence is determined by the flow 
$ t . This allows t to be included in the canonical relation A of the scattering operator Tr a F a , which is a 
map to spacetime distributions. Parameterized by y, rj and t, A becomes (51 ). The injectivity follows 
from the parameterization. 



Lemma 1. Let tq — T s (xo) and s(x) = tq — T s (x). // J(x, £) 
that maps (xo,£ ) onto itself. Its Jacobian is 



□ 



$ s (x, £) then J is a smooth bijection 



det(9( Xi€) J(x ,£ ) 



1 



l£ol 



n s (x ), 



(59) 



which is nonsingular by the DSE. 



Proof. For x in the neighborhood of Xq one has s(x) € /, so Q s is defined. The smoothness of J follows 
directly from the smoothness of x n- T s (x) and <& s in its arguments including s. The Jacobian results 
from the straight forward calculation 

<9(x,£) J ( x o>£o) = d (Xi £)$ (xo,£o) + 9 s $ o(x ,lo) ® d (Xi £)s(x ) 




c(Xo, 



l€ol 



^0 / 
' ICol 

|l |a xC (x ) 



-|^ol 9 xc(x ) / 

n s (x ) 
)9xT s (x ) 



(-5 x t s (x ) o; 



Herein we substitute the right-hand side of the characteristic ODE (11) for 9 S $ S 



□ 



4 Reverse time continuation from the boundary 

The receiver wave field is modeled by the reverse time continued wave u T . In this section, we show 
that u T is the result of a pseudodifferential operator of order zero acting on the continued scattered 
wave Uh- We refer to it as the revert operator P. 

The processes that are modeled by P are the propagation of the scattered wave field from a certain 
time, say t — t c , to the surface at x n — 0, the restriction of the wave field to the acquisition domain, 
the data processing, and eventually the continuation in reverse time. The revert operator suppresses 
the part of the scattered wave field that cannot be recovered because the contributing waves do not 
reach the acquisition domain. The data processing comprises a spatial smooth cutoff on the acquisition 
domain, the removal of direct source waves and the removal of waves reaching the surface following 
grazing rays. The final reconstruction represents a field related to bicharacteristics that intersect the 
acquisition domain M only once, and in the upgoing direction. 

Let u be the solution to the homogeneous wave equation. When we apply the result of this 
section to develope the inverse scattering, we set u = u n - Let M be a bounded open subset of 
{(x, t) 6 W l+1 | x n = 0} and let Tmu denote the restriction of u to M. We denote x' = {x\, . . . , x n -\), 
so (x',£) are coordinates on M. The field u r is an anticausal solution to 

[c(x)- 2 d t 2 - A] u r (x, t) = 8{x n )F M T M u{x! , t), (60) 
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where Fm is defined as follows. 

The boundary operator Fm consists of two types of factors. A pseudodifferential operator accounts 
for the fact that the boundary data for the backpropagation enters as a source and not as a boundary 
condition. This operator is given by 



- 2LD i c-yi - c 2 Di 2 Dl„ D t = P 1 ^, D x , = i" 1 ^. (61) 

The singularity in the square root is avoided by the cutoff for grazing rays, see below. The second 
type of factor is composed of three cutoffs: 

(i) The multiplication by a cutoff function that smoothly goes to zero near the boundary of the 
acquisition domain. The distance over which it goes from 1 to in practice depends on the 
wavelengths present in the data. 

(ii) The second cutoff is a pseudodifferential operator which removes waves that reach the surface 
along tangently incoming rays. Its symbol is zero around (x',f, £ ,w) such that 

c(y', 0)|(7/, 0)| =±w, 

and 1 some distance away from this set. If, given the velocity and the support of Sc, there are 
no tangent rays, this cutoff is not needed. 

(iii) The third cutoff suppresses direct rays. Since the velocity model is assumed to be known, these 
can be identified. 

We write £ ,w) for the symbol of the composition of these pseudodifferential cutoffs. The 

principal symbol of Fm is then 



- 2iwc- 1 V /l - C 2 W - 2 £' 2 f A f(x',t,£». (62) 

The decoupling procedure presented above yields two fields u a and Ub, associated respectively with 
the negative and positive frequencies in u. We will show that u?^ and u r ^ depend locally on u a and 
Ub in the following fashion, 

X u I , a (;t) = X [P a .(t)u a {:t) +R 1 {t)u h (;t)] and 

(63 

X«r,b(-, t) = X [Pb{t)u h (;t) + R 2 (t)u a (;t)] . 

Here, P a (t) and Pb(t) are pseudodifferential operators described below and R\{t) and f?2(i) are reg- 
ularizing operators, and x is a cutoff because the source in equation (60) causes waves in both sides 
of x n = 0. Note that the decoupling, which so far was mostly a technical procedure, turns out to be 
essential to characterize the reverse time continued field. The revert operator in matrix form will be 
defined as the i-family of pseudodifferential operators 

Waves are assumed to hit the set M coming from x n > 0. We assume supp(x) to be compact and 
contained in the set with x n > 0, and we invoke the following assumption 

bicharacteristics through M and supp(x) intersect M only once and with dx n /dt < 0. (65) 

The operators P a and depend on Fm and on the bicharacteristic flow in space-time between the 
hyperplanes t = and x n — 0. Let X s denote the set W 1 x {s} C K™ x K t . The bicharacteristic flow 
provides a map 

(X, ^ (y^x, €, t), t, T/' a (x, C, t), -C(x)|£|), 
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from T*X Q to T* M. The principal symbols of P a , Pb, which we will denote by p a , p\>, are then given 
by the following transported versions of ^ m'- 

U/ f Y ft / *M(yi(x,€,«),* > *I / a(x J €i*).-c(x)|€|) when 3< with y a (x,^,<) e M, 
*^,a(x,€) - | otherwise, (66) 

and ^x .h is defined similarly using the (yb^b) fl° w - We can now state and prove the following 



Theorem 3. Let it r . a; and u a , Ub, x an d M be as just defined. Equation (63) holds, in which P a 
and Pb are pseudodifferential operators in O'pS (M n ), the principal symbols of which are given by 

p a (*;x,£) = *x ,a(x,£) and pb(t;x,£) = *x ,b(x,C), (67) 



respectively. The operators P a , Pb satisfy property (4-9) as far as they are uniquely determined consid- 
ering the cutoff x in (63). 

The proof will be presented in the remainder of this section. If we take Cauchy values at t = t c , 
then for small \t — t c \, T]^u(x.',t) can be described by the local FIO representation of the solution 
operator. This representation can also be used for the description of the map from Tmu to «,■(•, i c ). 
The result can then be proven by an explicit use of the method of stationary phase. For longer times 
we apply a partition of unity in time to Tj^u(x! ,t), so that for each contribution the length of the 
time interval is small enough to apply the local FIO representation. Egorov's theorem will be used to 
reduce to the short time case. Alternatively one could consider one-way wave theory as a method of 
proof. 



Proof. We prove (63 1 for some given t. Without loss of generality we may assume that t = 0. The 
field u, by assumption, solves the homogeneous wave equation and is determined (possibly modulo a 
smooth contribution) by the Cauchy values ito, a = u a (-,0) and Mo,b = Wb(-,0). Consider the equation 
dtu a = ~\Bu a . In this proof we write S a (t, s) instead of S a (t — s) for the operator that maps initial 
values at time s to the values of the solution at time t. We write S a (-, s) for the operator that maps 
an initial value at time s to the solution as a function of (x, t), t > s. We will write S a (t, •) for the 
operator that gives the anticausal solution to (dt + iPu a )w a = f a , 



S a (t,-)fn = -J S a (t,s)f a (; S )ds. 

Note that S a (t, ■) maps a function of (x, t) to a function of x and that S a (t,-) — —S(-,t)*. The 
restriction operator Tm introduced above maps C°°(M. n xl)-) C°°(M) and is given by 

T M u{x',t) =u(x',0,t), (x',t)eM. 

The adjoint of this operator is given by the following. With auxiliary function / it is: 

T* M f{x,t)=5{x n )f{x!,t). 



These operators are well defined on suitable sets of distributions. We use the notation (cf. (60)) 

/ M (x',t) = F M T M u{x',t) 

and study the map (uq a , Uo.b) ^ fu- It follows from the results on decoupling that 

fid = F M T M {cS a u ^ + cS h u .b), (68) 

modulo a smooth error. Following this decoupling, we analyze the map u a t— > Pjw^M c 'S'a u o,a- 
To begin with, there exists a pseudodifferential operator Fm such that 

F M T M u = T M F M u 
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modulo a smooth function. This holds for a distribution u that satisfies |£| < C\uj\ in WF(«) for some 
C, like the solution of the homogeneous wave equation. Naturally, i*jvf(x, t, £, ui) ^ Fm (x', t, £ , w), 
because then the symbol property would not be satisfied around the line (£ ,u) = 0, £ n 7^ 0, but 
in the neighborhood of this line the symbol can be modified without affecting the singularities since 
|£| < C\uj\ in WF(u). Thus the first term in (68) is given by 

TmFmcS^-,0) (69) 

acting on uo.a, which is a product of Fourier integral operators. The operator S a (-,0) has canonical 
relation 

{((y a (x, £, t),t, r, a (x, £, i), -c(x)|£|), (x, £))}. (70) 

The operator i 7 ^/ removes singularities propagating on rays that are tangent or close to tangent to the 
plane x n — 0, and the restriction operator to x n = has canonical relation 

{((y',t,ri',u;),(y',0,t, V ',r) n ,oj))}. (71) 



As tangent rays are removed, the composition of canonical relations (71) and (70) is transversal. 
Therefore, (69) is a Fourier integral operator. Moreover, from assumption (65) it follows that the 
canonical relation is the graph of an invertible map, given by 

{((yl(x,^t),t ) r ? , a (x,^ ! t),-c(x)|^|),(x,0) | * s.t. y a (x,€,f) G M}, 

or more precisely a subset of this set, taking into account the essential support of Fm- 

Next, we consider the map /m >-> X u r,a- We insert a pseudodifferential cutoff S(x', t, D X ', F> t ). It 
cuts out tangent rays and is defined such that SFm = Fm- Using the decoupling procedure of section 



2.5 the source (/ a , /b) for the inhomogeneous wave equation are given by (/ a , /t>) T = A(0, c/) T , hence 



Xu r a satisfies 

X«r,a(- S 0) = X Sa(0, •) ( ^ c)T* M Ef M 

There exists an operator S such that T'^Ef = ET M f at least microlocally on the set |£| < C\oj\ for 
large C. Then XMr,a(">0) is given by the operator 

acting on /m, modulo a smoothing operator. 

The operator ^5^(0, -)(^B _1 c)S is a Fourier integral operator with canonical relation 

{((x,0,(y(x,^t),t,r ? (x,^t),-c(x)|^|))| |^(x,€,t)| >e, e > 0}. 



For an element (y', 0, rj', w) with \uj\ > c\r}'\ there are two rays associated, namely with rj n = ±-J 'c~ 2 w 2 — 
The + sign propagates into x n < for decreasing time, the — sign points into x n > 0. The contribu- 
tions are well separated because of the cutoff for tangent rays present in Fm- Because of assumption 
(65) and the cutoff x, the contributions with + sign can be ignored. We write si } (0, •)(-|S- 1 c)ST* / 
for the Fourier integral operator that propagates only the singularities from M into the x n > region 
for decreasing time. By a similar reasoning as above, this is a Fourier integral operator with canonical 
relation contained in 

{ ((x, £), (y'(x, £, t),t, r/(x, t), -c(x)|€|)) | y n (x, & t) = 0}. 

Again this is an invertible canonical relation. 

The next step is the composition of the maps (uo, aj u o.b) /m and /m > (w r ,a(*! 0), Ur,b( - i 0)). As 
both maps are Fourier integral operators with canonical relations that are the graph of an invertible 
map, the composition is a (sum of) well defined Fourier integral operators. The fields u a and u r a 
are associated with negative uj, Uh and u r b with positive ui. One can verify that the "cross terms" 
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M o,a >-> u r,b('jO) and M .b i-> u rj a(',0) are smoothing operators. The maps u ,a >-> u r.a(',0) and 
«o,b !-> Ur,b(-,0) are pscudodifferential operators. The principal symbol p a (0;x, is the product of 
tyx and another factor. 

We proceed under the assumption that \E'm(x / , t, t) is supported in the region < t < t\, with 
ii sufficiently small such that the explicit form of the Fourier integral operator can be used. This 
assumption will be lifted at the end of the proof. We treat only the map u , a >-> u r.a(',0), the map 
«o,b !-> Ur,b{-, 0) can be done in a similar way. The map wo, a >->• /m can then be written in the form 

/r,a(y',0) - J f a( fwd )(y',^x,Oe i(Q(y ^°^« ) - x -« ) U o,a(x)^dx, 

where the amplitude satisfies 

a( fwd V,t,x,£) = -2ix(ar n ) Wv /l - c(y', 0) 2 w" V 2 ^det^x) * M (y, t, r)', to) mod S a (R 2n x R n ), 

(72) 

where co — d t a = — c(x)|£|, 77 = <9 y a and det(<9 y x) is the Jacobian of the ray flow as explained earlier. 
The adjoint of the map Jm Xu r ,a(-,0) is given by E*TMC^B^ 1 S a (-,0)x, and is a Fourier integral 
operator with the same phase function a(y', 0, t, £) x • £ and amplitude 



a( bkd V,M,0 = ix(^)(-^- 1 )c(y) V /dct(a y x(y',t,C))S mod S- 2 (K 2 " x R n ). (73) 
The map /m >->• X u r,a(-,0) is therefore given by, with the notation z instead of x € R™, 

«r,a(z,0) - ^ HI a(bkd) (yT 4) z , C)e i(-a(y',0, t ,C)+-C) /Af (y ' ; t) dC dy > dt 

Therefore, the map u , a h-> x( z n)«r,a(-, 0) has distribution kernel K(z,x) given by 

1 1 f a o*<D(y', t, z, C) a (fwd) (y', t, x, $)c i (- Q (y'- - t -«+ Q (y'-°- t -«+-<—«d y ' dt ^ dC (74) 

Using a smooth cutoff the (£, £) integration domain can be divided into three parts, one with \C\ < 2|£|, 
one with |£| > and a third part containing (£, £) = (0,0). In the first part, the method of 

stationary phase can be applied to the integral over (y',£,0 using |£| as large parameter. We show 
that there is a function g(z, x, £) such that 

^bkd) a (fwd) e i(-a(y',0,t,C)+a(y',0,t,C)+-.C-x.O ft = fl ( Zj Xj £) e *(— *)* , ( 75 ) 



(2*)' 



and such that g(z, x, £) is a symbol that has an asymptotic series expansion with leading order term 
satisfying g(x, x, £) = (x, £)■ 

The first step in this computation is to determine the stationary points of the map 

$ : (y', t, C) -a(y', 0, t, C) + a(y', 0, i, £) + z • C - x • £. 

By the properties of a, g^y$ = if and only if g^y(y',0,C) = ^%y(y',0,|) if and only if 
(y',0, t, C) and (y',0, i, £) are associated with the same bicharacteristic and hence C — £• Requiring 
that the derivative with respect to £ is gives that 

-%*(y,t,£)+z = 0. 

Therefore, the bicharacteristic determined by (z, £) must be the same as the bicharacteristic determined 
by (y',t, £). Let ^(y'l C; x j £) be a C°° cutoff function that is one for a small neighborhood of 
(y')*iC) around the stationary value, and zero outside a slightly larger neighborhood. From the 
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lemma of nonstationary phase one can derive that the contribution to g from the region away from 
the stationary point set is in S~°° . 

At this point, observe that the second part, with \C, \ > can be treated similarly, with the role 
of C and £ interchanged. In this case the stationary point set is in the region where the amplitude 

x E"). The third part, C,£ 



is zero, and its contribution is of the form (75), but with g in S~°°(R : 
around zero, also yields such a contribution with g E S^ oa (M 2n x W 1 ). 

To treat the case (y',i, C) around the stationary point set, we apply a change of variables in the 
phase function. Setting y n = 0, it can be written as 



(y,*,€)-o(y,*,C)= J ^<x(y,tX + s(t-Q)ds = (t-Q- J ||(y,t,C + *(e-C))«k 



1 da 



The goal is to rewrite ( [74| by change of variables into 
1 



(27T) 



^ a (bkd) a (fwd) e i((€-c)-x+z-c 



dX 



d(y,t) 



dC dX, 



(76) 



so the next step is to prove that d ^ y , ^ is an invertible matrix at the stationary points. It is clear that, 
with £ = £ and y n = 0, 

da 

%'U,0^(y,U)-x( y ,«), 



where (y, t, £) >->• x(y, i, £) was discussed in subsection 2.4 The matrix |p is non-degenerate. Then 
we apply the implicit function theorem to the map x (y',t) obtained by setting y' = y'(x, t) in 
which t is such that y n (x,t) = 0, and use that there are no tangent rays, to obtain that the matrix 



dx 

d(y',t) 



has maximal rank at the stationary points, while the Jacobian satisfies 



dy 



dy n 



dt 



The integral ( 76 1 has a quadratic phase function £ • (X — z) , and can be performed as usual in the 
method of stationary phase [THl lemma 1.2.4] This shows that <?(z,x, £) satisfies the symbol property. 
Using ([72]) and (Pf3l) it follows that 



ff (x,x,0 = -2i(-c|£|)v/l-c 2 u;-V 2 



dx 
dy 



* M (-i) C ( 5 c)(c|€|)- 1 c(y) 



9x 
dy 



dx 
dy 



dyn 
dt 



Two terms need to be worked out, namely \J\ — c 2 to 2 r]' 2 = cos(#m), in which 9m is the angle of 
incidence of a ray at M, and \^§r\ — c(y) cos(^jvf)- Therefore indeed we have 



<?(x,x,£) = * Xo mod S X (I 



')• 



This concludes the proof of the small time result. 

Next we extend this to the result for longer times. By a partition of unity we can write ^5>m as a 
sum of terms with t € [s, s + ti] for some s. It is sufficient to prove the result for each term, and we 
may therefore assume t G [s, s + tx] in the support of ^ m- By a change of variable t to t — s, it follows 
that 

P a (s) d = S<->(«, •) H^-V) T; / F M T M 5 a (-, S ) g OpS°(]R") 
with principal symbol 

$ Xs (x,£) = $ M (y'(x,$,t - s),t,ri'(x,Z,t - s),-c(x)\Z\), with t s.t. y n (x,£,t - s) = 0. 



19 



From the group property of the S a (t, s) it follows that 

P a (0) = X Sa(0, S )P a (s)S a ( S ,0). 

The evolutions operators S a (0,s), S a (s,Q) are each others inverses. According to the Egorov theorem 
[43| section 8.1] the operator P a (0) is a pseudodifferential operator. For the symbol we find that it is 
given by (x, £) H> ^ Xs (y( x , s ): >?(x, s)), i.e. by \&x - This completes the proof. □ 



5 Inverse scattering 



This section deals with the inverse scattering problem. The diagram in figure [2] shows how we theoret- 
ically approach RTM. The forward modeling is given by r — > u — > d in the diagram. The reflectivity 
function r causes a scattered wave field u, giving the data d by restriction to the surface x n — (recall 
x' = [x] i :rl _i). The bottom line of the diagram shows the inverse modeling. Data d is propagated in 
reverse time to the reverse time continued wave field u r . This wave field is mapped by the imaging 
operator G to the image i. The resolution operator R is the map from the reflectivity to the image as 
result of the forward modeling and the inversion. The scattering operator F maps the reflectivity to 
the continued scattered wave Uh- As explained, this field can be seen as the receiver wave field in an 
idealized experiment. It contains all rays that are present in the scattered wave, regardless whether 
they can be reconstructed by RTM. The revert operator P removes parts that are not present in the 
receiver wave field. The field Uh, central to the analysis, is not actually computed. 



We obtain the main result, the imaging condition (89 1, in two steps. We propose the imaging 



operator G and show in theorem [5] and its proof that it is a FIO that maps the reverse time continued 
wave field to an image of the reflectivity. Hence it is an approximate inverse of the scattering opera- 
tor. From this operator we subsequently derive an imaging condition in terms of solutions of partial 
differential equations, g and u r . We first discuss a simplified case with constant coefficient. 



Instead of condition (651 we have the following condition for the RTM based inversion 



bicharacteristics that enter the region x n < do not return to the region x n > 0. (77) 
This will ensure that u r is properly defined for the purpose of linearized inversion. We also recall 



the assumption that there is no source wave field multipathing, formalized as the property (50). The 



assumption that there are no direct rays from the source to a receivers is incorporated in P, i.e. by 



means of ^m, cf. (62 1. 




u(x, t) 



i(x) it r (x, t) d(x', t) 



Figure 2: Diagram showing the theoretical approach to RTM. 



5.1 Constant background velocity 

In this subsection we consider the case of constant background velocity c with a planar incoming wave, 
propagating in the positive x^ direction. The scattered field will be described by 

[c- 2 d? - A] u(x, t) =A5(t- c- 1 x 3 )r(x), (78) 
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which is a slight simplification of ( 39 ) . For simplicity the analysis will be 3-dimensional, but it applies 
to other dimensions as well. 



The solution of the PDE (78 1 is given in the (£, t) domain by 



e M€|(*-s) _ e -ic|€|(t-a) 



2ic|C| 



f(£,s)ds, 



(79) 



where, for now, we denote by / the right hand side of (78). The Fourier transform of / is hence 
needed. Let r(£i, £2, £3) be the Fourier transform of r with respect to (sci,ara) but not £3. The Fourier 
transform of A5(t — — )r(x) is given by 



^A5(t-^)^ 1 ,&,x 3 )dx 3 = cAe- i & ct r(Z u &,ct). 



(80) 



Next we use (79) and (80) to solve (78), and we make a change of variable cs — z. This yields 

i-tc 



2ic|€| 



Ae-^ 2 r(Ci,^,z) dz. 



We can recognize in this formula a Fourier transformation with respect to z. However, the Fourier 
transform of r is not evaluated at £3, but at £ 3 ± |£|, because z occurs at several places in the complex 
exponents. Under the assumption that the support of r is contained in < X3 < ct (in other words, 
that we consider the field at time t such that the incoming wave front has completely passed the 
support of the reflectivity), the formula equals 



(81) 



The field in position coordinates is given by the inverse Fourier transform of this, i.e. by 
1 



u(x, t) 



(2tt) e 



M\ct^± ?{i + (0j 0j m _ e -m<**A?(£ _ (o, 0, |£|)) 



2ic|£| 



2ic|£| 



e ix ^ d£ (82) 



The two terms yield complex conjugate contributions after integration. To see this, change the inte- 
gration variables in the second term to — £, and use that the property that r(x) is real for all x is 
equivalent to r(£) = r(— £) for all £. Therefore 



u(x, t) = 



1 



(2tt) 3 



: icA 



Re / e- i l€l d+ «-«^ff(e-(0,0,|€|))de. 



(83) 



There are three wave vectors in (83), £ is the wave vector of the outgoing reflected wave, (0,0, |£|) 
can be interpreted as the wave vector of the incoming wave, while £ — (0, 0, |£|) can be interpreted as 
the reflectivity wavenumber, which, for a conormal singularity for example, would be normal to the 
reflector. 

In this simplified analysis we assume that the reverse time continued receiver field u r satisfies a 
homogeneous wave equation with equal final values (after the scattering) as u, like u\ x in (40), i.e. it 



results from an idealized experiment as explained in section |3.1| This means that u r is also given by 
(83), except that this formula is now valid for all t. 

The basic idea of imaging is to time-correlate the source field with the receiver field. Approximating 
the source field by AS(t — X3/6) this becomes evaluating the receiver field at the arrival time of the 
incoming wave and multiplication by A. Hence, a first guess for the image would be Iq = j4«(x, xs/c). 
This, however will not yield an inverse. Using some advance knowledge we will define instead as our 
image 

(84) 



2 

J(x) = -2-r(9t + cd X3 )u(x,x 3 /c). 
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We have from ^83 1 
2 



c2A (d t +cd X3 )u(x,t) {2n)3 



Re 



1 - 



e- i|€|et+bt ' € f(€-(o,o,|e|))d€. 



Setting t — x-ijc we find 
J(x) = 



(2tt): 



Re 



1 - 



3 ix.(«-(o 1 o,l€l)f^_( 0)0) |^|))^ 



We carry out a coordinate transformation, 

C = €-(o,o,|€|), 



9£ 



= l - 



(85) 



(86) 



(87) 



The image of this transformation is the halfplane £ 3 < 0, while the Jacobian is as given in (871, and 
exactly equals the factor 1 — -fe from the derivative operator dt + cd X3 . Therefore by a change of 

variables (86) equals , 2 \ 3 Re J? <Q e lx '^r(£) d£. This can be rewritten as 



7(x) 



(2tt) 3 



e ix| f(|)rfe 



The right-hand side is almost the inverse Fourier transform, except for the exclusion of the set £3=0 
from the integration domain. This expresses the difficulty with inverting from direct waves. This 
simple calculation gives the motivation for the imaging condition (89 1 below, in particular, for the 
term involving the gradient <9 x u r (x, u). 



5.2 Imaging condition 



We present the main result of the paper. The imaging condition yields a mapping of the source wave 
<?(x, t) and the reverse time continued wave u r (x, t) to an image i(x) of the reflectivity. We will show 
that the following imaging condition yields a partial inverse, 



i(x) = 



1 

2n 



iw|g(x,w)| 2 



g(x,w)u r (x,w) - 



c(x) 2 



<9 X ?( X 7 W ) • <9 x tl r (x,w) duj 



(89) 



Here fi(cj) is a smooth function, valued on a bounded neighborhood of the origin, and 1 outside a 
slightly larger neighborhood. These neighborhoods are obtained in the proof of the theorem. 

To characterize i(x), the relation (48) between £ and £ is important. We observe that the inverse 
function of (48) is defined on the halfspace 



{CeM"\0|C-n s (z) <0}. 



(90) 



The function p a (T s (z); z, £(C))> Pa the principal symbol of the revert operator, is in principle defined 
only on (90). However, due to the DSE, it is zero for £ near the boundary of this halfspace and we 
will consider it as a function on W l \0 that is zero outside (90). With this definition, the function 
(z, £) i-> p a (T s (z); z, £(£)) is an order symbol. 



Theorem 4. Let image i(x) be defined by (89), and assume (31), (38) and (77). Define operator R 
by the map from the reflectivity r to the image, i?r(x) = i(x). Then R is a pseudodifferential operator 
of order zero, and its principal symbol satisfies 



p.s.(ii)(z,C)=Pa(T 8 (z);z^(C))+p a (T s (z);z^(-C)) 1 
where the map (z,£) 1— > p a (T B (z);z, £(C)) * s as just described. 



(91) 
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Operator R will be referred to as the resolution operator. From the proof of the result it can be 



seen that the first contribution on the right-hand side of (|91|) corresponds to the negative frequencies 

the first, of these 
Hence, the map 



and the second contribution to the positive frequencies. As the supports, i.e. (90) for the first, of these 

< R n \0. 



two terms are disjoint, (91) defines a symbol that is one on a subset of . 



d i— > i given by (89) can rightfully be called a partial inverse. 

The imaging condition (89) is based on the actual source field g. Before proving theorem |4j we 
derive an intermediate result with an imaging condition based on the source wave traveltime T s (x), and 
the highest order contribution to the amplitude ^4 S ( X )- Let w € £'(YxR) be an auxiliary distribution. 
Let operators H and K be defined by 



Hw(y,t) 



1 



d t 2 [d t + c(y)n s (y) ■ d y ]w(y, t), 



(92) 



My) 
Kw(z)=w{z,T s (z)). 

Operator K is a, restriction to a hypersurface in K n+1 . Operator H is a pseudodifferential operator. 

—5+1 . ... _ _n±l 

Operator d t 2 is to be read as the pseudodifferential operator with symbol uj i-> d{ui)(iui) t~ in 
which a is a smooth function, valued 1 except for the origin where it is 0. Because P and F are defined 
as matrix operators, we define V\ — (l 0) which projects out the first component of a two-vector. 
We define the imaging operator G = KH . 



Theorem 5. If (31), (38) and (71) are satisfied and R is given by 



Rr = Gu T = GViPFr, 



(93) 



then R is a pseudodifferential operator of order zero with principal symbol given by (91) 



Proof. We first work out the details for the negative frequencies, leading to a characterization of 
i? a = KHcP a F a . We then consider the positive, and add the contributions, R = R a + R\,. 

(i) We show that the composition R a = KHcP a F a is a FIO and that it is microlocal, i.e. has 
canonical relation that is a subset of the identity. The kernel of operator K is an oscillatory integral, 



Kw(z) = (2tt)" 



-y)+MT.(.)-t) w ( y; t ) d ( Vj t ) dT] duJ 



(94) 



with canonical relation 



T = {((z,0),(y,t,ff,w)) | (y,r?)ery\0, t = T s (y), ,el\{0}, 

z = y, 6 = q + c(y)- 1 u J n s (y)}. (95) 

First consider Kir a F a , which is the composition of K and ir a F a with canonical relations given re- 



spectively by T (95) and A (51). We consider the composition of the Fourier integrals K and 7r a F a , 
using the composition theorem based on the canonical relations, see JTHl Theorem 2.4.1] or [35]. Let 
((z,9), (x, £)) e T o A then there exist a (y,Tj) € T*Y\0 that is not in V Sj t, time t = T s (y) and 
uj = — c(y)|?7| such that ((z, 6), (y, t, r\, ui)) € T and ((y, t, 77, ui), (x, £)) € A. As a result one has 
(x, ^) = < l > T B (x)-T s (y) (y ; I 7)i which means that x and y are on the same ray and separated in time 



by T s (y) — T s (x). Condition (50) (SME) now implies that this ray must coincide with a source ray. As 
source rays are excluded, i.e. (y, rj) ^ V s t, the only possibility is that x = y. The conclusion is that 
(z,0) = (x,C). 

It is straightforward to establish that the composition of canonical relations is transversal, and that 
the additional conditions of the composition theorem of FIOs are satisfied. Hence Kir a F a is a FIO with 
canonical relation contained in the identity. The operators H and P a are pseudodifferential operators, 
and 7r a and P a can be constructed such that WF(P a w) C WF(7r a u>) for all w. The conclusion is that 
R a = K HcP a F a is a FIO with identity canonical relation, and hence a pseudodifferential operator. 
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(ii) We show that R = KHViPF is a pseudo-differential operator that can be written as the integral 



(1091 below. For F we use the local expressions (42). Because P is a t- family of pscudodiffcrcntial 



operators and F a p is a i-family of FIOs, the composition P a (t)F a (t)p is a FIO with phase inherited from 
F a (t)p, i.e. ip T . The highest order contribution to its amplitude is p a (t; y, 9y(/3 T )Ap. The composition 
with H can be done similarly, because P a (t)F a (t)p can also be viewed as a FIOs with output variables 
(y, t). In this proof we will denote the highest order contribution to the amplitude of H cP a (t) F a (t) p 
by Ahpf(v, ti, x, £). It can be written in the form 

a / , t\ {-i , c (yK(y) ■ 9 y a \ 2ic(y)p a (f i; y, d y a)aA s {x)d t a 

A HPFiy,ti,x,€) = I 1 H . (9bJ 

V d t a J c(x)A s (y) 

For all occurences of a and a the arguments are (y, t\ — T s (x); £). 

Next we consider the applicaton of the restriction operator K. We have already argued that R a is 
a FIO with canonical relation contai ned in the identity. This implies that, to prove the theorem, it is 



sufficient to do a local analysis using (42 1. The local analysis shows again that R a is a pseudodifferential 
operator, but also gives the required explicit formula for the amplitude. 

The local phase function of KHcP a (t)F a (t)p will be denoted by ip(z,x,£). Applying K to ip T , i.e. 
setting t = T s (z), yields 

V>(z, x, £) = a(z, T s (z) - T„(x); £) - £ • x. (97) 
The stationary point set of ip, denoted by 'J, is given by the triplets (z,x, £) that solve 

%*(z,T s (z)-T B (x);0=x. (98) 

The interpretation of (z, x, £) G $ is that a ray with initial condition (x, £) arrives at z after time 
lapse T s (z) — T s (x). Application of the SME and the DSE now implies that z = x. 

Below we will define a transformation of covariables. To prepare for this, we introduce a smooth 
cutoff function x '■ Z x X x M. n \ {0} — > K accordingly. A Fourier integral may be restricted to a 
neighborhood of the stationary point set at the expense of a regularizing operator. Therefore, x( z , x , 
is set to 1 in the neighborhood of \t and elsewhere. This means that x is close to z in supp(x). 
The second issue is related to the DSE, which is required for the definition of the transformation. The 
cutoff x is assumed to also remove singularities on a neighborhood of the direct rays. We set x( z i x > £ ) 
to if £ lies within a narrow conic set with solid angle f2(z) around the principal direction n s (z). The 
solid angle f2(z) will be discussed later. We can hence write 

R a r(z) = (2tt)- b J J e^^ x '«)x(z,x,|)A H PF(z,T s (z),x,0^r(x)dx. (99) 

in which, of course, the integration domain is implicitly restricted to supp(x). 

Next we introduce covariable 9 to transform phase rp into the form 9 ■ (z — x). By definition 
9(z, x, £) = — J Q d x ip(z, x(/x), £)d(x in which x(/i) = z + /i(x — z). The phase function now transforms 
into 



Jo 



(100) 



To better understand the transformation and to determine the new domain of integration, i.e. <?(supp(x)), 
and the Jacobian we apply the chain rule to the definition of %p. This leads to 

0(z,x,€)=€+ / 9 t a(z,T s (z)-T s (x);|)5 x T s (x)d M . 
Jo 

There exists an x such that (z, x, £) 6 Tt s (x)-t s (x)i 1 - e - x an d z are connected by a ray. Note that x = 
xr(z,T s (z) — T s (x),£) will do, see section 2.4 for notation xp. By using the identities d t a = — c(x)|£| 
and c(x)9 x T s (x) = n s (x), one gets 

e(z,x,€)=€-|€|n(z,x,€) with n(z,x,0= f'^-n^dp. (101) 
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The Jacobian now follows from this result. By an easily verified calculation, one finds 



\d € e\ 



det I 



n(z,x,£) 



n(z,x,£) 



(102) 



With these formulae at hand a sensible choice can be made for the solid angle f2(z). The angle 
must be large enough to meet the following inequality for all elements of supp(x): 



|£-n(z,x,£)| < |£|min{l,|n(z,x,0| 2 }. 



(103) 



We will now give the motivation. For £ n- 0(z,x, £) to be injective, given (z, x), the Jacobian must 
be nonzero. This is true due to the inequality, which is nontrivial if |n| > 1. This affirms the local 
invertibility, and an easy exercise proofs its injectivity. A second argument concerns the domain of 
integration 0(supp(x)). The inequality guarantees that 0(z,x, £) • n(z, x, £) < for all points in 
supp(x), which is nontrivial if |n| < 1. This fact will play a role in gluing i? a r and i?br together, 
which will be done in following paragraphs. Because x is in the neighborhood of z, so are x and x. 
This implies that n(z,x, £) is close to n s (z), and |n| sa 1. This is as close as needed by narrowing the 
spatial part of the cutoff function x around the diagonal of ZxX. 

By using the new variable i? a r's integral expression (991 transforms into 



R a r(z) = (2tt)- 



(supp(x)) 
-1, 



A 6 (z,x,0)e ie(z - x) d0r(x)dx, 



(104) 



where wc define 

Ag(z, x, 9) = \d € e\- l X (z, x, C) A HPF (z, T s (z), x, £)• (105) 

Concerning the integration domain it can be observed that, for a given (z,x) the set 0(supp(x)) is 
contained in the halfspace {9 € M™\{0} | 9 ■ n< 0}. 

(iv) While the expression ( |104 ) defines a pseudodifferential operator of order 0, it is given in a 
non-standard from. It differs from a regular pseudodifferential operator, because the the amplitude 
As(z, x, 9) depends on (z, x, 9) and not only on (z, 9). Another amplitude that does not depend on x 
can be found by 



n -i 

A^(z,z,0) + V / D eh d Xk A^(z,z + //(x - z), 0) d/z, 



(106) 



which is an application of formulae (4.8)-(4.10) of Tre ves 45J. The first term is the principal symbol 
of i? a , which has order 0. The second term in (106) does not contribute to the principal part, it 
corresponds to a pseudodifferential operator of order — 1. We will denote by As(z, 0) (with two 
arguments) the symbol of R. 

To evaluate of Ahpf (96) on the diagonal one applies ffify, the relation d t a(z,Q;£) = — c(z)|£| for 



the phase and the result a(z,0;£) = 



2c(z)|C| 



for the amplitude. This yields 



A H pf(z,T s (z),z,£) = ( 1 - ) Pa (T s (z);z,|) = \d € 9\ Pa (T s (z); z, £) 



(107) 



see also (102). In view of (105 )-( 1071, we have p.s.(Ag)(z, 0) = x(z,z,£)p a (T s (z);z,£). Note that 
r] = d y a £ holds on the diagonal, and that £ = £(0). 

We now come back to the formal role of cutoff function x- By requiring x( z j z jO = 1 on 
supp(p a (T s (z); z, £)) the cutoff function can be left out. This requirement is allowed because Q(z) 
in the construction of x can be chosen arbitrarily tight by narrowing the spatial support of x around 
the diagonal. Therefore 

p.s.(A g )(z, 9) = p a (T s (z); z, £(0)). (108) 

(v) A key step is the inclusion of both negative and positive frequencies. In section [2] we saw that 
<x(x, t; £) e lA "( x '' ; ^ and 6(x, t; — e lA ^( x ' t;- £) have a symmetry relation: They yield complex conjugate 
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contributions (note the — signs). The consequences of this property can be traced through this proof. 
We find that Rt,r{z) — (27r)~™ //_ e ( supp ( x )) Ag_(z, — 9)e l0 '^~^d9 r(x) dx, and consequently, modulo a 
regularizing contribution, 



Rr(z) = (2tt)" 



A ft (z, 0) + A g (z, -0)1 c i6 {z - x) d8 r(x) dx. 



(109) 



The ©-integration is over the full space because the definition ofAs (z, 0) can be smoothly extended 
such that it is zero outside the domain 0(supp(x)). In view of (108) this proves the claim. □ 



Proof of theorem^ The first step in deriving the imaging condition is to rewrite operators H , K and 
G (92). Let w(y:,t) again be an auxiliary distribution. In this section u>(x,w) will denote its temporal 
Fourier transform. Because u>(x, t) = J e lut w(x, uj) dui, one has 



Kw(x) = 



1 

2n 



[iuj + c(x)n s (x) • d x }w(x, w) 



(110) 



Applied to the reverse time continued wave field u r (x, uS), equation (93) becomes 



1 
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AuT, 



^ a, \ 0^) [iw + c(x)nix) ■ <9 x ]u r (x, w) dw. 
A s (xJ 



(111) 



The next step is to eliminate T s (x), A s (x) and n s (x) by expressing them in terms of the source 



field explicitly. The principal term of the geometrical optics approximation of the source ( 34 ) is 



Function cr, introduced in (34), is smooth and has value 1 except for a small neighborhood of the origin 



where it is 0. Later we will examine the effect of the subprincipal terms of the source and the division 
by its amplitude. One naively derives the following identities 



3 iwT 8 (x) 



1 



A,(x) 



cr(w) 



c(x)n s (x) 



(iw) 2 ?(x,w) 

c(x) 2 d x ff(x,o;) _ c(x) 2 d x g(x, u) 
— iwg(x, w) 



(112) 



i<x>g(x, w) 

in which it is used that the second equation is real-valued. Substitution of involved factors occurring 



in the integral (111) yields 



i?r(x) 



1 

2^ 



g(o;)cr(a;) i c(x) 2 d x g(x, w) • <9 X 
kjg(x,w) 



u r (x, w) dw. 



(113) 



(iw) 2 g(x,w) 

We will finally argue that the division by the source amplitude is well-defined and that the sub- 



principal terms in the expansion for <?(x, u) do not affect the expression for the principal symbol (91 ). 
The source wave field is free of caustics by assumption. The transport equation yields that, on a 
compact domain in spacetime, there exists a lower bound L > for the principal amplitude, thus 
\Aq (x, x s , oj)\ > L. Division by A is therefore well-defined, and from its homogeneity and the inequal- 

j4(x,x s ,ui) 



ity (33) it can be deduced that there exists a constant C > such that 



A (x,x a 



1 



< 



C 

i+M' 



For 



|w| sufficiently large, division by A is therefore well-defined. We choose 1—0 wide enough such that 
all u € supp(J7) are high and satisfy d{uS)a{uS) — 1. The difference between ^- and jjj is of lower order 



in u. By construction it holds that Aq(x,x s ,oj) = A s (x) (icu) 2 on supp(f2). Taking (113) we replace 
da with Q, to define the imaging condition (89). □ 
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6 Numerical examples 



In this section, we give numerical examples to support our theorems. The general setup of the ex- 
amples was as follows. First a model was chosen, consisting of a background medium c, a medium 
perturbation (contrast) Sc — cr, a domain of interest and a computational domain. The latter was 
larger than the domain of interest and included absorbing boundaries. Data were generated by solv- 
ing the inhomogeneous wave equation with velocity c + Sc, and a Ricker wavelet source signature at 
position x s = (0, 0), using an order (2,4) finite difference scheme [IS]. The direct wave was eliminated. 
The operator (61) could be applied in the Fourier domain since in the examples c was constant at 
the surface. The backpropagated field was then computed using the finite difference method, and the 



same for the source field. Finally the imaging condition (89) was applied to obtain an approximate 
reconstruction of Sc. 

As we mentioned, only a partial reconstruction of Sc is possible in realistic situations. Relation 
(48) and the wave propagation restrict the directions of £ where inversion is possible. The frequency 
range present in the data also restricts the length of £, according to (48) and using that |£| = c _1 |w|. 



To be able to compare the original and reconstructed reflectivity we used bandlimitcd functions for 
Sc, which where obtained by multiplying a plane wave with a window function. Such functions are 
localized in position, by the support of the window, and in wave vector by the plane wave. 

Our first example concerns a gradient type medium with c(x\, X2) = 2.0 + O.OOIX2 with c in km/s 
and X2 in meters. Our model region was the square with x\ and Xi between and 2000 meters. 
The purpose was to show a successful reconstruction of velocity perturbations at different positions 
and with different orientations in the model. We therefore chose for Sc a linear combination of three 
wave packets at different locations, with central wave vector well within in the inversion aperture. We 
included one with large dip, as one of the interesting abilities of RTM is imaging of large dips. The 
results of the above procedure are shown in Figures [3] and [4j The reconstruction of the phase is 
excellent. However, the reconstructed amplitude is around 8-10 % smaller than the original amplitude. 
Possible explanations for this are inaccuracies related to the linearization and to a limited aperture. 



dc (km/ s) 



(a) 




reconstructed dc (km/s) 



(b) 



Figure 3: Example 1: Velocity perturbation and reconstructed velocity perturbation. The background 
medium is a gradient c = 2.0 + 0.001x2, with x 2 m meters and c in km/s. 

Our second example concerns a bandlimited continuous reflector. For a continuous reflector one 
might expect less loss in amplitude when compared to the localized velocity perturbations. One of the 
strengths of RTM and wave equation migration in general is that multipathing is easily incorporated, 
where in our case of single source RTM, multipathing is only allowed between the reflector and the 
receiver point. To see this in an example we included in our background model a low velocity lens 
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dc and reconstructed dc at xl=400 m 



dc and reconstructed dc at xl = 1400 m 



dc and I'ccon^Tructcd dc at xM=600 ir_ 




Figure 4: Example 1: Comparison of some traces from Figure [3] at x\ = 400 m, x\ = 1400 m and 
X2 = 600 m. 



at (800, 1200) m. The background medium including some rays, as well as some data are plotted in 
Figure [5] The velocity perturbation was located at x 2 = 1600 m. The results of this example are 
given in Figure [6] The reconstruction of the phase is again excellent. The amplitude varies somewhat 
depending on location, being about 0-10 % too low. The smooth tapering which was applied has 
diminished smiles and amplitude variations, but not fully eliminated them. The multipathing leads 
to singularities in the inverse of the source field u^., around (£1,2:2) = (1900, 1000) m, which leads to 
the two artifacts that can be seen there. 

Velocity, rays and fronts 




x (m) w sim data, direct arrival removed 

(a) (b) 



Figure 5: Example 2: (a) A velocity model with some rays; (b) Simulated data, with direct arrival 
removed. 



7 Discussion 

We presented a comprehensive analysis of RTM-based imaging, and introduced an imaging condition 
condition involving only local (data point and image point) operators which yields a parametrix for 
the single scattering problem for a given point source. 

We make the following observations concerning our inverse scattering procedure: (i) The symbol 
of the normal operator associated with a single point source contains a singularity which has been 
observed in the form of "low-frequency" artifacts [IHIISDI 19, 48, 23]. Our imaging condition yields a 



parametrix and naturally avoids this singularity, (ii) The square-root operator (61), a factor of Fm 
introduced in section|4j can be removed with dual sensor (streamer) data, that is, if the surface-normal 
derivative of the wave field is measured. We note that Fm is available only microlocally. (iii) Division 
by the source field, in frequency, can lead to poor results when its amplitude is small. There are two 
main reasons why this can occur. First, a realistic source signature can yield very small values for 
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dc (km/s) 02 reconstructed dc (km/s) 

(a) (b) 



Figure 6: Example 2: (a) Velocity perturbation; (b) Partial reconstruction of the velocity perturbation. 

particular frequencies in its amplitude spectrum. Moving averaging in frequency typically resolves 
this situation [531 IH1 03] ■ Secondly, the illumination due to propagation in a velocity model of high 
complexity may result in small values; spatial averaging over small neighborhoods of the image points 
may be benificial. (The cross-correlation imaging has been adapted by normalization with the source 
wave field energy at the imaging points as a proxy to inverse scattering [HE].) 

The acquisition aperture, and associated illumination, is intimately connected to the resolution 
operator R. This operator is pseudodifferential and the support of its symbol expresses which parts of 
the contrast or reflectivity can be recovered from the available data. Partial reconstruction is optimally 
formulated in terms of curvelets or wave packets. A detailed procedure, making use of the fact that 
the single scattering or imaging operator is associated with a canonical graph, can be found in [15] : 
see also [15] , 

We have addressed the single-source acquisition geometry, which arises naturally in RTM. One can 
anticipate an immediate extension of our reconstruction to multi-source data, but a major challenge 
arises because the single source reconstructions are only partial. Because each of the single source 
images result in reconstructions at different sets of points and orientations, in general, which are 
not identified within the RTM algorithm, averaging must be excluded. However, techniques from 
microlocal analysis can be invoked to properly exploit the discrete multi-source acquisition geometry. 
(We note that in the case of open sets of sources the generation of source caustics will be allowed.) 
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